Mass transport subject to time-dependent flow with non-uniform sorption in porous 

media 



OS 

o 
o 

(N 



Marie- Christine Neel/'Q Andrea Zoia, 2 'Q and Maminirina Joelson 1 '^ 

1 Universite d 'Avignon et des Pays de Vaucluse, UMR 1114 EMMAH, 84018 Avignon Cedex, France 
2 CEA/Saclay, DEN/DM2S/SFME/LSET, Bat. 454, 91191 Gif-sur-Yvette Cedex, France 

We address the description of solutes flow with trapping processes in porous media. Starting from 
a small-scale model for tracer particles trajectories, we derive the corresponding governing equations 
for the concentration of the mobile and immobile phases. We show that this formulation is fairly 
general and can easily take into account non-constant coefficients and in particular space-dependent 
sorption rates. The transport equations are solved numerically and a comparison with Monte Carlo 
particle-tracking simulations of spatial contaminant profiles and breakthrough curves is proposed, 
so to illustrate the obtained results. 
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I. INTRODUCTION 



Contaminant migration in porous media is often char- 
acterized by non-Fickian (anomalous) transport: fol- 
lowing injection, the spread of the pollutants plume 
might grow nonlinearly in time, (x 2 (t) — (x(t)) 2 ) ~ 
i 7 , 7 7^ 1, and the resulting concentration profiles 
display a non-Gaussian behavior [H, H, [I]- In con- 
trast, particles flow in perfectly homogeneous media 
(where the Fickian advection-dispersion mechanisms ap- 
ply) gives rise to linear spread and Gaussian shapes: see, 
e.g., [J] and references therein. Many concurrent pro- 
cesses may explain the observed deviations from Gaus- 
sianity. For instance, the presence of irregularities at 
multiple space scales 0, Q, the complex structures of 
flow streams [7|, [8| and saturation/stagnation distribu- 
tion within the medium Q, and the physical-chemical 
or bio-physical exchanges of the pollutant particles with 
the surrounding material [lol [Tl| make the homogeneity 
hypothesis questionable. 

Anomalous transport often displays non-universal fea- 
tures: different physical conditions lead to concentra- 
tion profiles that, while sharing some properties (such 
as the scaling law for the spread, for instance), can 
not be interpreted within a single coherent framework. 
This is especially apparent in presence of boundaries [l2[ , 
and explains the coexistence of several models aimed at 
understanding and predicting solutes dynamics in com- 
plex materials. Among them, some of the most widely 
adopted formulations are the Continuous Time Random 
Walk (CTRW) [IL [land the fractal Mobile-Immobile 
Model (f-MIM) HElli: both have indeed been applied 
with success to the analysis of experimental data ranging 
from laboratory to field scale [J 0, [TH, [l6| ■ In particu- 
lar, these approaches are well suited to shed light on the 
'heavy tailed' (power-law decaying) breakthrough curves 
(BTCs) that are frequently measured at the outlet of ex- 
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perimental setups, and the non-Gaussian shapes of spa- 
tial contaminant profiles. For a detailed discussion on 
the distinct features, adva ntag es and limitations of these 
models see, e.g., [1 [lS [3 HI E3 • 

The long-time asymptotic behavior of these formula- 
tions may look very similar (up to an appropriate renam- 
ing of the coefficients) |l) when the analysis is limited to 
given physical quantities (such as BTCs), whereas rel- 
evant discrepancies may appear at a closer inspection 
(by examining, e.g., spatial profiles as we will see in 
the following). Moreover, the governing equations for 
anomalous transport sometimes appear in the literature 
under different forms, equivalent when all parameters, 
e.g., the dispersion coefficient or the velocity, are con- 
stant and uniform. This equivalence may break down 
when the parameters depend on position and/or time 
(see, e.g., Pi|). In all such cases, one must adopt more 
precise hypotheses on the microscopic solute dynamics 
in the traversed medium, so to single out an appropriate 
model for the experimental data under consideration. In 
general, however, solutes trajectories are hardly accessi- 
ble by experiments (at least in the context of contam- 
inant transport in porous media), so that one must re- 
sort to ensemble-averaged macroscopic measurable quan- 
tities in order to discriminate between hypotheses. For 
instance, the comparison between model-predicted and 
experimentally measured BTCs and spatial profiles may 
allow choosing the most appropriate conceptual frame- 
work. 

According to the f-MIM model, in the hydrodynamic 
limit the evolution of the contaminant plume density 
is ruled by a transport equation involving integral- 
differential operators of non-integer (fractional) order in 
time, as shown in [l7|, 13 on the basis of numerical 
and theoretical arguments. Here the hydrodynamic limit 
refers to the fact that we are observing the plume be- 
havior at space and time scales much longer than those 
characterizing the typical particles displacements. The 
prototype equation for the evolution of so-called frac- 
tional dynamics is the Fractional Fokker-Planck Equa- 
tion (FFPE) [H SJ: indeed, f-MIM and FFPE share 
many features, and show a similar asymptotic behavior. 
However, these two equations are not equivalent, and rep- 
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resent hydro-dynamic limits of distinct small-scale models 
for particles trajectories. 

In particular, the FFPE corresponds to random walks 
performing Gaussian jumps (in potential fields) that take 
random durations to be completed [2l|, |22j, |23[ , whereas 
f-MIM describes random walks involving immobile peri- 
ods (of random duration) and Gaussian displacements at 
each time step [l4], EH, E3, H3] . Conceptually, the former 
is an expedient means of describing a broad spectrum 
of velocities, such those characterizing flows in hetero- 
geneous and/or nonsaturated media, whereas the latter 
allows distinguishing between a solid matrix (where par- 
ticles are stuck, such as in a low permeability region: the 
so-called immobile phase) and a bulk flow (where par- 
ticles undergo advection and dispersion processes: the 
so-called mobile phase). 

The aim of our work is to provide a generalization of 
the f-MIM approach to the case of non-constant flows, 
and space-dependent sorption rates, which can com- 
monly arise in transport experiments. Non-uniformity 
can occur in both field and laboratory scale measures, 
and may involve also sharp discontinuities in the physical 
properties of the traversed media. Based on a small-scale 
description of particles trajectories, we will derive the 
corresponding governing equations for the macroscopic 
quantities, namely the mobile and immobile densities, 
in such a way that non-constant coefficients are easily 
taken into account. In the case of constant velocity field 
and uniform sorption, similar results were obtained via 
subordination theory [ijj], along the strategy proposed 
in [22|, HE HH, [2lJ • In the case of non-constant coeffi- 
cients, stronger arguments have to be used. Finally, in 
order to corroborate the proposed results, we will com- 
pare Monte Carlo particle-tracking simulations of solutes 
spatial profiles and BTCs with the numerical solutions of 
the governing equations. 

This paper is organized as follows: in Sec. [II] we recall 
the small-scale model for flow with trapping processes 
in porous media, on the basis of the f-MIM formalism, 
and provide an extension to non-constant flows. Then, in 
Sec. Mil we derive the densities of the mobile and immo- 
bile contaminant phases at small scale, and in Sec. lIVI the 
corresponding governing equations in the hydrodynamic 
limit. These equations are discretized and solved nu- 
merically, and the obtained solutions are compared with 
Monte Carlo simulations results in Sec. [V] Conclusions 
are finally drawn in Sec. IVII 



II. A SMALL-SCALE MODEL FOR FLOW 
WITH TRAPPING PROCESSES 

As customary, we begin by representing the stochastic 
trajectory of a contaminant particle in a porous medium 
as a random walk x £ t ' T undergoing advection and disper- 
sion. Superscripts £, r denote the characteristic length 
and time scales, respectively, of the process. First, we 
briefly recall the essential features of the standard Gaus- 



sian models that usually describe small-scale displace- 
ments of a contaminant plume in homogeneous saturated 
materials. Then, we focus on heterogeneous and/or un- 
saturated media, where the broad distribution of per- 
meabilities and different flow regions experienced by the 
tracers is mirrored in the possibility of trapping events 
at each visited spatial site, the walkers dynamics being 
otherwise similar to that observed in homogeneous mate- 
rials. These trapping events affect the sojourn times and 
thus alter the typical scales of average displacement (i.e., 
velocity) and spread (i.e., dispersion) of the contaminant 
plume. 



A. A Gaussian model for homogeneous materials 

For homogeneous saturated materials, it is usually pos- 
sible to identify an average flow field v(t), so that at each 
time step [t — r, t] particles are advected over a distance 
fj,(t) = f v(t')dt' . Dispersion is usually taken into ac- 
count by adding random (symmetrical) jumps of char- 
acteristic length scale £ to the advective contribution /i. 
The total displacement during r can be therefore written 
as Ax = n + ££, where £ is a random number drawn from 
a probability density function (pdf) y>i(£) with zero mean 
and unit variance. It follows that ipi(-) — £~ 1 (pi(-/£) is 
the density of ££. Usually, one assumes that <pe(£) is a 
normal pdf with zero mean and standard deviation equal 
to £, which means that the typical scale of fluctuations 
around the average particle displacement is £. It is well 
known that the scaling (hydrodynamic) limit of such (in- 
dependent) random walks is attained when r — > in such 
a way that the ratio £ 2 /2t converges to some limit D. 
The parameter D defines the dispersion coefficient. Un- 
der this scaling limit and when v is constant, the stochas- 
tic process x e t ' T asymptotically approaches the Brownian 
Motion (BM) x t with drift, whose concentration P(x,t) 
(i.e., the probability density of finding a walker at posi- 
tion x at a given time t) is shown to satisfy the following 
Fokker-Planck Equation (FPE) [H[ 

8 t P(x, t) = d 2 x DP(x, t) - d x vP(x, t). (1) 

Note that the FPE (fTJ) is equivalent to the Advection- 
Dispersion Equation (ADE) when D is uniform [2^ |. 
When D depends on x, Eq. (fTJ) still represents the hy- 
drodynamic limit of random walks as above, satisfying 
I 2 1 It = D(x) [2!|. For sake of simplicity, in the fol- 
lowing we will refer to the case where D is uniform, i.e., 
D(x) — D. By virtue of the Central Limit Theorem, the 
results recalled above actually apply more generally to a 
broad class of random walks where tpi(£) is a generic sym- 
metric jump length pdf with finite second moment. Af- 
ter a sufficient number of displacements, these processes 
all converge to BM in the hydrodynamic limit (provided 
that the variance of the process is equal to £ 2 , and with 
constant v). The fact that BM is the basin of attraction 
of a large spectrum of random walks, independently of 
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the specific choice of the jump length pdf, can explain 
the success of Eq. ([T]) in interpreting experimental con- 
taminant transport data, at least limited to homogeneous 
saturated materials Q • Remark also that Eq. |T]) can be 
derived by building upon a constitutive relationship for 
the particles flux (probability current) !F e,T (x,t) of the 
process x l t ,T , as recalled in Appendix [Pi In the hydrody- 
namic limit, the flux converges to 

T{x, t) = vP(x, t) - d x DP(x, t), (2) 

which is the well-known Fick's law for homogeneous 
D. Then, mass conservation principle d t P(x,t) — 
—d x J 7 (x,t) yields Eq. ([1]). This holds for infinite do- 
mains, or for domains limited by absorbing boundary 
conditions (walkers are removed upon touching the bar- 
riers) . For detailed accounts concerning the derivation of 
Eq. (J), see, e.g., [11 SHI- 

B. Dispersion with trapping events 

Suppose now that the traversed material is affected by 
small-scale heterogeneities. For instance, we might con- 
sider unsaturated porous media with variable saturation, 
where particles may be retained by stagnation regions [9( . 
We mi ght also think of biophysical effects in porous ma- 
terials [11], or chemical reactions taking place on the sur- 
face of a duct traversed by fluid flow. The homogeneous 
random walk proposed above is evidently inadequate to 
address such situations. From the point of view of mi- 
croscopic trajectories, a natural means of accounting for 
sorption is that walkers are given the possibility of being 
trapped at the end of each displacement, i.e., upon reach- 
ing a new spatial site. The trapping probability h{x) is 
for sake of generality space-dependent, since traps may 
be non-uniformly distributed. Previously, [ItI fl9L |22| ] 
considered constant trapping probabilities. We further 
assume that the sojourn time at the traps is itself a ran- 
dom variable t w . Experimental evidences suggest that 
these retention times lack a characteristic scale (i.e., their 
average is not defined), so that it is commonly assumed 
that t w obeys a power- law decaying pdf d, 01 • 

We introduce the scaled variable t w = t^'^W, where 
W is a random number obeying a pdf ip and 7 a scaling 
exponent, and denote by \& the associated survival 
probability, *(t) = ijj(t')dt' . The quantity *(f) 

expresses the probability that the trapping time is 
longer than t. It follows that the rescaled pdf of t w is 
■0 7 = r _1 / 7 -0(t/r 1 / 7 ), and the rescaled survival prob- 
ability is $7 = ^{t/r 1 / 1 ). We make now the following 
hypothesis for the pdf of the retention times: 

Hi) the pdf "0 is concentrated on R + , with survival 
probability of the kind *(t) = Xt^ /T(l - 7) + K(t), K 
being a function integrable over R + , with < 7 < 1. 

Assumption Hi is satisfied by pdfs whose asymptotic 
behavior is a power-law. For instance, we might consider 



Pareto laws [3l|, |32j , or maximally skewed Levy laws with 
exponent 7, which are concentrated on R + precisely for 
0<7<lpll3ll.l35l l36l|. Intuitively, the exponent 
7 quantifies the degree of heterogeneity of the porous 
media: small values of 7 denote strong deviations from 
the usual Gaussian transport model, i.e., anomalously 
long retention times. 

While the retention times correspond to the immobi- 
lization of the walkers (immobile phase), a more precise 
description of the time spent during displacements (mo- 
bile phase) is needed. Several scenarios can be conceived, 
depending on the time of occurrence of the dispersive 
jump within a mobile period [t — T,t\. For instance, the 
jump could take place instantaneously at the beginning 
of the period, or at the end; or it could occur at a ran- 
dom time, uniformly distributed in [t — T,t]. Also, we 
could imagine that the jump is not instantaneous, and 
takes the whole time span [t — r, t] to be completed. On 
the other hand, the endpoints of successive displacements 
do not depend on the considered scenario, nor are trap- 
ping events affected. As shown in Appendix [3] all these 
possible small-scale random walks converge to the same 
scaling limit when t, £ — > 0. Then, for convenience we will 
focus on the simplest case: we assume that walkers per- 
form a single instantaneous dispersive jump during the 
time interval r, taking place at the end of each mobile 
period. 

Finally, for sake of generality, we also introduce a (pos- 
sibly time- and space-dependent) source term r(x, t), rep- 
resenting tracer injection. 

In the following, we will show that in the hydrody- 
namic limit the walkers density P(x, t) for the process 
described above satisfies 

8 t P = %DU x „, h P - d x vHx, lth P + r. (3) 

In Eq. ([3]), the non-local in time operator TL\ n ^h is 
the inverse of the (also non-local in time) mapping 
Id + \h(x)I ~+ , which entails the fractional integral of 

order 1 — 7, namely 7 whose definition is recalled 
in Appendix [B] Here Id denotes the identity opera- 
tor and A > is a constant parameter. In fact, as 
shown in [li^ ]. Tix^ji is the time convolution of the ker- 
nel ^Ei_ 1 [-Xh(x)t 1 ~'J}, where E a is the Mittag-Leffler 
function described in [33, Hg, [39[ . In [li| it was shown 
that Eq. ([3]) governs the evolution of the particles con- 
centration for constant v and uniform D, with h(x) = 1, 
building upon the results of for this case, and as- 
suming r(x, t) — 5(t)6(x), Eq. ^ is equivalent to the 
fractal MIM model 

(d t + \dl) P{x, t) = d x {Dd x - v)P(x, t) (4) 

introduced in [3], where is the Caputo derivative of 
order 7 (see Appendix [Bj . The exponent 7 < 1 char- 
acterizes the asymptotic behavior of the trapping times 
pdf (in Hi), and also the scaling of the plume spread; in 
this sense, 7 is the signature of the anomalous transport 
process. Indeed, the solutions of Eq. ((3]) have been shown 
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to decrease at large times as t 7 [14| , which could possi- 
bly explain the long tails experimentally observed, e.g., 
by and 0. In unbounded domains, and with constant 
and uniform coefficients, the spatial moments of the so- 
lute concentration in Eq. ((4]) were shown to decrease as 
powers of time related to the exponent 7 [171 ]. 

As a special limit, setting 7 = 1 in Eq. (0| yields the 
well-known MIM model with retardation factor 1 + A [24j . 
which corresponds to tracers experiencing random re- 
tention periods with finite characteristic (mean) dura- 
tion, comparable to the time spent in the average flow 
field. Nevertheless, the solutions of the MIM model or 
the ADE |T]) fall off much more rapidly than any power 
of t and are thus inadequate to interpret experimental 
data showing heavy tails such as those of [zL H|. We 
will see further below that, rather than representing the 
hydrodynamic limit of random walks satisfying H\ with 
7=1, the MIM model corresponds to a pdf ip with a 
finite average. 

The parameter A determines the relevance of the reten- 
tion mechanism with respect to the Fickian transport: 
when A = 0, all equations above collapse to the ADE. 
Moreover, it provides the scaling parameter for ip, and 
carries dimensions of a power 7 of time. This is easily 
seen in the Laplace space, where ip(s) ~ 1 — As 7 in the 
limit s — *■ 0, according to H\. 

III. PROBABILITY DENSITIES FOR THE 
MOBILE AND IMMOBILE PHASES AT SMALL 
SCALE 

Particles performing such random walks with sorption 
can conceptually be separated into two distinct 'phases': 
at each time step, walkers that are trapped are said to 
be in the immobile phase, whereas walkers that are not 
are said to be in the mobile phase. In the following, 
we proceed to derive an explicit relation that links the 
particles densities in the two phases, for definite values 
of length- and time-scales i and r. The hydrodynamic 
limit will be addressed in next Section. 

Let Pf' T (x,t) be the density of trapped particles, at 
location x at time t, and P^ l T {x 1 t) the density of mobile 
walkers. In order to establish the desired relation be- 
tween the two spatial densities P-' T and P^i T , we make 
use of the ancillary pdfs pj' T of just arriving at point x 
at time t, and p£ of just being released by a sorbing site 
at time t. 

Except just after having been injected into the system, 
mobile particles at position x at time t have two alterna- 
tives. Either they may have completed a mobile period at 
time t — t' (0 < t' < t), without being trapped; or, they 
may have been trapped and then released, at a distance 
f. , v(8)d9 from x. Both possible events are followed by 
a convective displacement that may not be completed at 
time t. The displacement completed at time t has ampli- 
tude L(t,t') = f*_ t , v(6)d0. Remark that L(t,t') = vt' if 
v is constant. Hence, we have the following probability 



balance 

P^ T {x,t)= [ T Tt>y L(t ,t>)[f' T + r}(x,t)dt', (5) 
Jo 

where the quantity 

f t > r (x,t)= J fc(x,t) + [1 - h(x)]pf T (x,t) (6) 

is the pdf of just beginning a mobile period at time t and 
position x, after a previous mobile period (i.e., particles 
just injected by the source are excluded). Convective 
displacements are represented by means of the operators 
T u and y w , which denote translation in time and space, 
respectively; i.e., T u G(t) = [HG] (t — u), H being the 
Heaviside step function, and y w g{x) = g(x — w). Eq. §5§ 
corresponds to scenario (SI) of Appendix [Al the disper- 
sive jumps occurring at the end of each mobile period. 

Immobile particles that are in x at time t must have 
jumped there previously, been trapped and stayed there 
up to t. Hence, denoting time convolutions of functions 
in R+ by *, i.e., F * G(t) = f* F(t - t')G(t')dt', we have 

P^(x 1 t) = h(x)* 1 *pf T . (7) 

To complete the mass balance above we need another 
equation. Particles just arriving at x at time i > r may 
i) have jumped at the previous time step without being 
trapped, ii) have been trapped and released, or Hi) may 
have been injected into the system by the source, in each 
case at time t — t. Hence, for pj' T (x,t), which appears 
on the right-hand side of Eq. , we have 

pf r (x,t)=T T y HttT) [f t 'T + r]*<pt (8) 

for t > t, where ★ denotes space convolution, i.e., 
/ * g(x) = J R f(x — x')g(x')dx' . Note that the explicit 
dependence of the convolution product on the variables 
(x, t) has been omitted. 

According to Eq. (0, r^ 1 P^ T (x,t) is the average of 
1~t'yL(t,t') [f e ' T + i"]( x ,t) over an interval of amplitude r, 
hence approximates < ^-3^i(t lT )[/ £ ' T + i"](x,t) when r be- 
comes small, at least for smooth functions of time. Since 
convolutions as in Eqs. ([7]) and ([S]) have a smoothing 
effect, this latter assumption is not necessary to ensure 
that replacing T T ~y L ( tiT )[f l > T + r) by t~ 1 P^ t results into 

a small error for P/' T . To check this argument, let us just 
split ^-3^L(t,T)[/ £:T + r ] into T~ 1 P J ^ r and the remainder, 
and estimate the influence of this latter in p^ T (x, i) and 
Eq. 0. We obtain 

P? < T (x, t) = h(x) [R e ' T P^(x, t) + E l ' T (x, t)] . (9) 

Further below, we will address the limit of the operator 
R e ' T , and check that E e ' T tends to zero when £, t — > 0. 
We have set 

R l ' T g(x, t) = t- 1 [% * g * tp t ] (x, t), (10) 
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FIG. 1: Breakthrough curve !F(x,t) at the right outlet of a 
domain with length L = 1. The time-dependent velocity is 
v(t) = 2sin(i/t), with D = 1, h = 1, A = 1 and 7 = 0.5. 
A point-source is located in xo = L/2 at time t = 0. The 
left outlet has reflective boundary conditions, while the right 
outlet has absorbing boundary conditions. Symbols represent 
Monte Carlo simulation (dots v = 5, squares v = 25, and 
crosses v — 50), solid line numerical integration. 

E l ' T (x,t) = ^ 1 *tpi*£ e ' T (x,t), (11) 

and 

e^(x,t) = T T y L{ttT) [f^ + r](x,t) ~ = 
Jo {%y L {t,r) -Te T y L{t ,er))[f' T + r](x,t)M. (12) 

The equations above provide the link between the mobile 
and immobile walkers densities, at small scale. 



IV. GOVERNING EQUATIONS 

In this Section, we derive the hydrodynamic limit of 
the small scale processes described above, and illustrate 
the macroscopic governing equations. The starting point 
is the limit of Eq. ©. 



A. Mobile and immobile densities 

Our aim is to show that the mapping R e ' T , defined by 
Eq. (|10p . converges to a fractional integral, whereas in 
Appendix [Cl we prove that the quantity E e,T can be ne- 
glected, provided that P^ T and f l ' T converge sufficiently 
smoothly when £, t — ► 0. Indeed, the mapping R l T com- 
bines convolutions in time and space, with kernels t -1 \I/ 7 
and tpe. The space convolution with kernel y>£ converges 
to the Id operator [4(|. The time convolution with ker- 
nel r _1 ^E , 7 splits into the sum of a singular term, and a 



FIG. 2: Total concentration profiles P(x,t) in domain with 
length L = 1, at fixed times. The time-dependent velocity 
is v(t) = 2sin(^), with D = 1, h = 1, A = 1, v = 50 and 
7 = 0.5. Initial data and boundary conditions are as in Fig.[T] 
Symbols represent Monte Carlo simulation (squares t = 0.1, 
crosses t = 0.2, and dots t = 0.3), solid lines numerical inte- 
gration. 

mapping that vanishes when r — > 0. Indeed, hypothesis 
Hi imposes 

= \t-i/T{l - 7 ) + T- 1 K(t/T 1 ft), (13) 

K being an integrable kernel. Moreover, in view of Ap- 
pendix [Bl the time convolution of kernel \t r /r(l — 
7) is precisely the fractional integral A/g^ 7 . For 
the second term, we have ||T _1 iir(</T 1 / 7 )|| i i(^ + ) = 
t 1 /"<~ 1 \\K \\l 1 {r+), hence Young's inequality (see Ap- 
pendix [Cj implies that the convolution of kernel 
T^Ki-fr 1 ^) is a mapping of L x ([0, T],X) that vanishes 
when t — ► for < 7 < 1. Therefore, recollecting the 
previous results, in the hydrodynamic limit we have 

Pi{x,t) = Xh(x)I^P m (x,t), (14) 

and 

P m (x,t)=H x , 1 , h P(x > t). (15) 

Expressions (fT4")) and (fTS")) provide the governing equa- 
tions for the mobile and immobile particles densities, re- 
spectively, at the macroscopic scale. When the pdf ip 
satisfies ip ~ i _7_1 with 7 > 1, so that it has a finite 
average, the survival probability \& is integrable and we 
have A = / + °° tip(t)dt = / + °° #(i)di. Then, the time 
convolution of the kernel r _1 \E , (t/r) approximates the 
identity operator Id when £ — > 0. Hence, the scaling 
t w = tW leads to Pi(x,t) — Xh(x)P m (x,t). Moreover, 
Eq. (Ull) still holds with Hx, 1>hh = 1/[1 + Xh(x)]. Thus, 
when the sticking times have a finite average, we recover 
the standard MIM model [24|, with a retardation factor 
1 + Xh (provided that h is uniform). 
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FIG. 3: Total tracers concentration P(x, i) at fixed times, on 
a domain of length L = 1. In [0, L/2] we have h(x) — 0, and 
h(x) — 1 in [L/2, L]. A point-source is located in xo = L/4 at 
time t = 0. Both ends of the domain have absorbing boundary 
conditions. The simulation parameters are: 7 = 0.8, A = 1, 
v = 0.5, and D — 0.2. Curves are plotted at times t = 0.25 
(squares), t — 0.5 (crosses), and t — 0.75 (dots). 

B. Particles fluxes 

In Appendix [D] we show that the probability current 
T l ' T (x,t) can be written as v(t)P^{ T (x, t) + J r jj T (x, t), up 
to an additive contribution that vanishes when £, r — > 0, 
and J 7 ^ ) T (x,t) — > —d x DP m . Then, using Eq. (fl"5|) and 
the definition P = P m + Pi yields the explicit expression 
for the total tracers flux 

F(x, t) = vHx.^hP - d x DH x , 1<h P, (16) 

which generalizes Eq. {2} to spatially distributed trap- 
ping events and variable velocity fields: this expression 
actually represents the Fick's law, as applied to P m . 
Combining this equation with mass conservation finally 
gives Eq. ©. 

The consistency of Eqs. (fH]) and (fTB|) will be verified 
by showing that solutions to Eq. (J3j) indeed describe the 
density of a plume of walkers performing the random 
walks in Sec. UH and that the associated particles fluxes 
satisfy Eq. (JIBJ). 

V. NUMERICAL SIMULATIONS AND 
COMPARISONS 

In a previous work, some of the authors discussed the 
use of numerical schemes discretizing Eq. ([3]) for the case 
of constant advection field v, and unit probability of un- 
dergoing a trapping event at the end of each displace- 
ment, i.e., h — 1 [19]. In this particular case, Eq. ([3]) 
is equivalent to the widely adopted Eq. Q, whose Ca- 
puto derivatives can be discretized according to various 



FIG. 4: Mobile tracers concentration P m (x,t) on a domain 
of length L = 1, at fixed times, with h(x), initial data and 
boundary conditions are as in Fig. [3] The simulation param- 
eters are: 7 = 0.8, A = 1, v = 0.5, and D = 0.2. Curves 
are plotted at times t — 0.25 (squares), t — 0.5 (crosses), and 
t = 0.75 (dots). 

existing numerical schemes [U S3- An alternative inte- 
gration method was proposed in |19j |. so to take advan- 
tage of the conservative form of Eq. ([3]). This scheme 
can be easily extended to the more general situation ad- 
dressed here, i.e., the fractal MIM equation (J3J with time- 
varying velocity and spatially-dependent sorption prob- 
ability. Therefore, we proceed now to display numeri- 
cal solutions of Eq. ([3j), and to compare them to Monte 
Carlo simulations of the microscopic-scale random walks 
described in Sec. [ill Indeed, in the hydrodynamic limit, 
the mobile fraction of an ensemble of random walkers un- 
dergoing the stochastic process described in Sec. [TT| ap- 
proximates the quantity P mi whereas the immobile frac- 
tion approximates Pi. After briefly revising the essential 
features of numerical integrations and random walk sim- 
ulations, we will present comparisons, so to illustrate the 
theoretical results of Sec.|lVl i.e., the Eqs. (fT4"]) and (TT5"]) 
and the subsequent Eqs. ([3J and (flT)]) . In particular, we 
will focus on cases where (although D is uniform) the 
fractal MIM formulation (J3J is not equivalent to 

[d t + X%] P{x, t) = d x (8 x D - v)P + Hx n , h r (17) 

which is the version of the more popular fractional differ- 
ential equation (|3]), suitable to deal with general source 
rates r(x,t). Comparisons between pde and random 
walks were presented in [l7| for infinite domains; here, 
we focus exclusively on bounded domains. 

A. Numerical methods 

Numerical integration of Eq. ([3J can be based on an 
implicit method with centered finite differences schemes 
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for space derivatives, described in [lj|, when P(x,t) is 
smooth, which is the case if h(x) does not show discon- 
tinuities. The non-local in time mapping TLx.-yji is ap- 
proximated by inverting a discrete version of the integral 
operator [Id + XH(x)Iq ~^ 7 ] [Hj]. Fluxes of tracers are fi- 
nally given by applying Eq. (|16l) . 

The Monte Carlo particle-tracking approach to the 
fractal MIM model described above consists in com- 
puting the trajectories of a (large) number N of inde- 
pendent particles performing successive displacements, 
whose rules are defined in Sec. [Til More precisely, let 
us denote by x n the location after the n th displacement 
of a particle that originated in xq at time to. This walker 
leaves x n at time t n , and we have 

rt n +r 

Xn+l = X n + J v(t')dt' + V2Dt£, 

where £ is a random Gaussian number with zero mean 
and unit variance, and either 

tn+l =t n +T + T^W, 

with probability h{x n +i), or 

t n +l = t n + T, (18) 

with probability 1 — h(x n +±). 

For the case of constant v and uniform h = 1, the re- 
sults in [l9j show that random walk simulations are in ex- 
cellent agreement with numerical integrations of Eq. 
The same holds for fluxes as in Eq. (fTBj) . obtained from 
both methods. In the following, we focus our attention on 
time-dependent velocities, and non-uniform probabilities 
h(x), in one-dimensional domains [a^av] with various 
boundary conditions. 

B. Time-dependent velocity 

We perform comparisons for a periodical velocity v — 
sini/t and h = 1. In Fig. [JJwe display the breakthrough 
curve (i.e., the outgoing flux from the domain) as a func- 
tion of time. In the context of solutes transport in porous 
media, this physical quantity is the most easily accessible 
by experiments, either at the outlet of laboratory-scale 
column setups, or at boreholes (wells) for large field-scale 
measurements. We consider a Id bounded domain of 
length L = 1 , with a reflective boundary condition at the 
left, X\ = 0, and an absorbing boundary condition at the 
right, x r = L. Velocity is positive pointing towards the 
right, so that we measure the outgoing flux at the right 
outlet of the domain. A point source is set at the center of 
the domain, xq = L/2. An excellent agreement is found 
between Monte Carlo simulation results and numerically 
integrated equations. 

Then, to further substantiate this analysis, in Fig. [5] we 
display the spatial concentration profiles P(x,t) for the 
total tracers concentration (at fixed times). In physical 




FIG. 5: Breakthrough curve J-(x, t) at the right outlet of 
a domain with length L = 1. A point-source is located in 
xo = L/4 at time t = 0. Both ends of the domain have 
absorbing boundary conditions. The simulation parameters 
are: A = 5, v = 2, and D = 1. Symbols represent Monte 
Carlo simulations (7 = 0.8 squares, 7 = 0.5 crosses, 7 = 0.3 
dots), solid lines the corresponding numerical integration. 

terms, these curves allow quantifying the average dis- 
placement and the spread of an initially close plume of 
injected solutes. Again, a very good agreement is found 
between Monte Carlo simulations and numerically inte- 
grated equations. 



C. Non-uniform trapping probability 

In the context of underground contaminant migration, 
a space-dependent sorption probability h(x) may be ex- 
pedient to represent the transition between zones of low 
and high permeability, as well as an alternation of satu- 
rated and stagnant regions. Actually, a non-uniform dis- 
tribution of the trapping sites is expected to be the most 
common situation in geological formations and complex 
soils. Conceptually, the simplest case is given by an 
abrupt variation between h — and h = 1 in two adja- 
cent portions of a given domain: this would correspond to 
the traversed medium being homogeneous and saturated 
in the former region, and unsaturated and/or heteroge- 
neous in the latter, where retention and sticking effects 
dominate. Then, the plume migration is Fickian where 
h = (memory effects are absent, and transport is ruled 
by the standard ADE dynamics, as seen from Eq. Q) 
and anomalous where h = 1 (transport is ruled by the 
memory kernel contained in the fractional integral). In 
the following, we illustrate this case by considering a Id 
domain of length L = 1, whose left portion [0, Xd] is char- 
acterized by h = 0, and whose right portion [xd, L] is 
characterized by h = 1. We assume Xd — L/2 and set 
absorbing boundary conditions at both ends of the do- 
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main. A point-source is located in xq = L/A at time 
t = 0, i.e., the particles are injected in the homogeneous 
and saturated region. In Fig. [3] we display the spatial 
concentration profiles for the total tracers concentration 
P(x,t) (at fixed times). The most striking feature of this 
transport process is the appearance of discontinuities in 
the resident concentration profiles at the interface be- 
tween the two portions of the domain, whereas the pro- 
files displayed in [l9| for h = 1 are continuous. This 
behavior has already been reported elsewhere in the con- 
text of the CTRW formulation (see, e.g., the discussions 
in [H, IH, m, Hil), and may be understood in terms of 
the two layers having distinct apparent porosities. In 
fact, particles experience different trapping times in each 
layer, and the abrupt variation of the sojourn times dis- 
tribution at the interface ultimately gives rise to sharp 
mass accumulation at the interface, as the particles flow 
is somehow hindered when going from the region with 
short retention times to the region with long retention 
times 46]. 

While the implementation of random walks for this 
case is straightforward, some care is necessary in dis- 
cretizing densities for numerical integration. More pre- 
cisely, the total concentration P{x, t) is discontinuous at 
Xd, hence directly discretizing Eq. ^ is not convenient. 
Using instead Eq. (|15p . which links P and P m , and then 
Fick's law applied to P m , is much more expedient, be- 
cause P m is smooth. Combining Eqs. (fT5j) and (J3j) yields 

d t P m (x, t) = H^. h [-d x vP m + d 2 x DP m + 

r(x,t) - Xh(x) T ^ ) P m (x,0+)] (19) 

which is easily discretized following the same lines as 
in 0. 

In Fig. 3] we display the spatial concentration profiles 
P m (x,t) for the mobile tracers concentration (at fixed 
times). In this case, the curves are smooth across the in- 
terface. This is because the particles flux satisfies Fick's 
law, as applied to P m , and the ADE does not allow for 
concentration discontinuities at the interface. Indeed, the 
flux contrasts local variations of P rn . Suppose that the 
profiles of P and P m have slopes of different signs (e.g., 
positive for P m and negative for P) at x: then, the parti- 
cles flux through x is negative, and flattens out the spatial 
increase of P m . On the contrary, the spatial decrease of P 
has no direct effect on the flux, which therefore does not 
act on this quantity. For both total and mobile concen- 
trations, a very good agreement is found between Monte 
Carlo simulations and numerically integrated equations. 
Finally, Fig. [5] shows the breakthrough curves at the col- 
umn right outlet, for different values of the exponent 7 
in the absorbing region: again, good agreement is found 
between Monte Carlo simulations and numerically inte- 
grated equations. The asymptotic behavior described 
by [H| m infinite domains with h(x) = 1 is recovered: 
even with h(x) 7^ in some intervals only, P(x, t) ~ i -7 , 
P m (x,t) ~ £ _7_1 and ^(Xjt) ~ t^ 1 ^ 1 when t — > +00. 
While density profiles obtained for a given value of 7 



with arbitrary h(x) show qualitative differences, fluxes 
look similar. Hence, BTCs alone are not enough to dis- 
criminate between cases. 



VI. CONCLUSIONS 

In this work, we have discussed a model of contami- 
nant particles flow with trapping events in porous me- 
dia. Building on the framework of the fractal MIM 
model, which describes advection-dispersion processes 
with sticking events in homogeneous flows, we have 
considered time-varying velocities and space-dependent 
sticking probabilities. We have first derived the small 
scale particles dynamics, on the base of a functional 
relationship between the densities of trapped and non- 
trapped walkers. This relationship stems directly from 
the asymptotic behavior of the trapping times distribu- 
tion, and gives rise to a modified Fick's law with memory 
for particles fluxes. Then, recalling mass conservation 
principle, we have obtained the corresponding governing 
equations for the evolution of the mobile and immobile 
phases densities. 

These equations have been derived by considering the 
hydrodynamic (scaling) limit of the underlying micro- 
scopic stochastic processes, i.e., by letting the space and 
time scale of the particles displacements be vanishing 
small, while preserving the macroscopic dispersion and 
advection coefficients. The transport equations, which 
contain non-local in time kernels in the form of fractional 
integrals, have been discretized and solved numerically 
by resorting to ad hoc algorithms. Finally, in order to 
corroborate our results, the contaminant concentration 
profiles and the breakthrough curves thus obtained by 
numerical integration have been compared with Monte 
Carlo particle-tracking simulations. 

The relevance and broad applicability of the trans- 
port equation ([3]) for the case of non-constant param- 
eters have been emphasized in both theoretical deriva- 
tions and numerical examples. In particular, we have 
addressed the case of time- varying velocity fields v(t) 
and space-dependent trapping probabilities h(x). In fact, 
the method developed here is more general, and may ap- 
ply also when the equation coefficients (e.g., h) depend 
on the densities of trapped and mobile walkers, which 
would result in a nonlinear version of Eq. ([3]), similarly 
as in [13, 13. 

Further extensions of our work will address the coexis- 
tence of several kinds of traps within the same porous 
medium, each trap being characterized by a distinct 
sticking time pdf. This approach would then give rise to 
slightly more complex mappings Tix.-yji with fractional 
integrals of distributed order |49fl . The simplest case 
would correspond to two kinds of traps occurring with 
probabilities hi and h,2, respectively, and sticking time 
pdfs -0 7l and ip~ 2 (satisfying hypothesis Hi with 71 < 72). 
Such a model could represent, e.g., multiple phases or 
regions, with distinct retention properties. Then, the 
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residence times of the solutes would be governed by the 
mapping T~i\ 1 ,\ 2 , 11 ,-y 2 ,hiM> defined as being the inverse of 
Id + A 1 /i 1 (2;)/q^ 71 +X 2 h2(x)lQ T 72 . Moreover, while at in- 
termediate times the solute dynamics is rather involved, 
at late times the total density would asymptotically de- 
crease as t~ 71 , i.e., transport would be dominated by the 
'slower' retention process. 



APPENDIX A: SOME SCENARIOS FOR 
DISPERSIVE JUMPS 

During a given mobile displacement between times t—r 
and t, we have a single dispersive jump, whose length is a 
random variable obeying if£. At the scale of microscopic 
particles trajectories, diverse scenarios may be conceived, 
which we denote by label (Si). Let xf' l ' T \uj) be the 
associated walkers paths. 

The source of randomness in each sample uj arises from 
the series of successive dispersive jump lengths J n {n> 1) 
and from that of trapping times T n ; for convenience, we 
set T n = if there is no trapping event after the n-th mo- 
bile period. Concerning dispersion, we might consider 
instantaneous jumps, occurring at the end of each mo- 
bile period (51), at the beginning (52), or at a random 
time uniformly distributed in the interval [t — r, t] (53). 
Alternatively, the dispersive jump might be thought of 
as being distributed along the total displacement, taking 
therefore a time [t — r, t] to be completed (54). 

Each u> corresponds to two sequences of numbers, that 
are the values drawn for the J n and T n . Of course, 
these points are identical for all the trajectories xf^' T%> (oj) 
started from xq at time t = 0. All trajectories pass 
through points (t n + r,x n +i) and (t n+ i,x n+ i), with 
t n +i = t„ + r + T n and x n+1 = x n + J n + J t " T v(9)d9: 
trapping periods correspond to segments beginning at 
point (t n +r, x n +i), and ending at (t n +T+T n , x n +i), that 
are common to all scenarios. Hence, the immobile walk- 
ers density Pf' T (x,t) does not depend on the scenario. 
Moreover, \x[ 1 ' 1 (u>) — {ui)\ = when t belongs 

to a trapping period, and \x[ l ' e ' T \uj) — {ui)\ < | J„| 

when t belongs to the n-th mobile period. Jumps J„ obey 
tpe, and (weakly) converge to zero when t — > 0. Hence, 
theorem 3.1 of [5(| implies that, if one among the pos- 
sible processes x[ l ' e ' T ^ converges to the scaling limit x t 
when i — > 0, then the same holds for all the other paths 
x t 

Hence, it follows that the walkers density P(x, t) does 
not depend on the specific scenario in the hydrody- 
namic limit, nor does the immobile walkers concentration 
Pi(x,t). Then, the mobile P m and immobile densities Pi 
do not depend on the scenario. 



APPENDIX B: FRACTIONAL INTEGRALS AND 
DERIVATIVES 

The fractional integral Iq + J of order a > is 

'§:+/(*) = fir l\t-t'Y- l f{t')dt\ 

r (c>0 Jo 

which generalizes the usual multiple integrals to non- 
integer order [4(| [El|. Observe that Iq + is bounded in 
LP[0,T] for 1 <p < oo HI. 

The Caputo fractional derivative d"f of order n < a < 
n+1 appearing in Eq. ^ is defined by 

n being an integer [13, [H, H(| • 



APPENDIX C: ESTIMATES AND LIMITS 

In this Appendix, we will prove that E e ' T — > when 
1,t — > 0, provided that P m and j t - T converge. We will 
make use of some technical results, which will also be 
used later in Appendix [Dl for fluxes. 



1. Hypotheses 

We will need some regularity assumptions for the 
source rate r, the velocity v and the densities P^ r and 

Hypothesis H 2 ) r(x,t) is the time derivative of some 
function p{x, t) S L 1 (R + , X), i.e., r(x, t) = dtp(x, t), and 
v is uniformly continuous. 

Observe that initial data of the kind P(x, 0+) = Po(x) 
lie within this assumption, with r(x,t) — 8(t)Po(x). 

Hypothesis H 3 . i) when £, r — * with D = £ 2 /2t, 
the density P^ T converges to P m in L 1 ([0,T],A), 
and d x P m belongs to this space, ii) The distribution 
f e ' T is the time derivative of some F l ' T that belongs 
to L x ([0,r],X): f> T = d t F e ' T : and F^ T — > F in 
L 1 ([0,T],X) when 1,t -> 0, with / = d t F. Moreover, 
Hi) d x {F z ' r + p) tends to d x (F + p) in L l ([0,T],X). 
Observe that point Hi) is not needed if v is constant. 



2. Statements 

We will make use of Youn g's inequality, which we 
reproduce here for convenience [52j: 
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Young's inequality. Let 1 < q',q,q" < oo, with X/q' + 
l/q" = 1 + 1/ q. Then, for F G L«'[Q,T] and G G L 9 "(i?), 
we have 

\\F*G \\l<i(R) < Il-Pll L?'[0,T] II ^lUv'CR)) 

and for F G £ 9 '(i?) and G G £ 9 "([0, T], Q) we have 

ll-F* G|| L ,([ , T],Q) < H- F llL<!'flll G llL'j"([0,T],Q)' 

where Q is a Banach space. 

We will prove the following proposition: 

Proposition 1. Suppose that hypotheses Hi, H% 
and i?3 are satisfied. Then, i) \& 7 * (pg * e^ 1 " — * in 
^([O,! 7 ],^). Moreover, ii) [H$(j)\ *e e ' T -> in the set 
<S'([0, T], X) of tempered distributions. 

Note that point z) implies i? £,T — > in the hydrody- 
namic limit. The proof will use the following lemmas. 

Lemma 1. i) Let w be a continuous function. 
Then, for g in i 1 ([0, T],X), Tt>y w (t')9 is a continuous 
function of if, with values in i 1 ([0, T], X). Moreover, 
m) the mapping Tt'y w {t') is a contraction in i 1 ([0, T], X). 

Lemma 2. If g G L x ([0,T],X), then 

Jo [^L(t,r - T eT y Htter) ]g{x,t)d9 in ^(M,*). 
when t — > 0. 

Consequence. If hypotheses #2 and i?3 are satisfied, 
then ti%y L{t , T ) ~ %ry L (t,9r)][F e ' r + P ](x,t)d0 -> in 
^([O,! 1 ],^) when 0. 



Moreover, since ?/» 7 is normalized, the Consequence and 
Young's inequality imply that the first term on the of 
r. h. s. of Eq. (|C2|) vanishes, which proves point i). 

For non-constant v, hypothesis H3 Hi) implies that 

$(■/£) *h = ^*/o [(«(*) - vit - T))T T y L(ttT) - 

(v(t) - v(t - 9T))T Tg y L(t , Te) ] [f' T + p] (x, t)d0 -> 

in i 1 ([0,T],X), since Ii — * while y>£ is normal- 
ized. Besides, W = * ft f£ (T T y L {t,T) ~ 
%eyL(t, T e)){F LT + p)(x,t)d9 may not belong to 
L x ([0, T],X). Nevertheless, U is the time derivative of 
$(•/£) * /^(T^i^) - r Te y L(tire) )[F^ T + p](x,t)d0, 
which vanishes in this space. This proves point ii). 

We have now to prove the Lemmas. In Lemma 1, point 
i) follows from the proof of Theorem 9.5 of [53| . stating 
that y a G is a continuous function of a, with values in 
L 1 (R), provided that we have G E L 1 (R). Point ii) is 
obvious. 

For Lemma 2, the function Tg T yL( t g T )g(x,t) be- 
longs to L 1 ([0,T],X), and depends continuously on 9 
by Lemma 1 i). Hence, it is Bochner-integrable (52| 
from [0,1] to Z^flp,! 1 ]^). Moreover, by Lemma 
1 m), T~Q T y L ( t Q T }g(x,t) — > g(x,t) pointwise, whereas 

l|r flT y L(tiflT)ff (i,i) - 5 ||f < Iblk m y = L^o.r],*) 

norm, so that dominated convergence proves the Lemma. 

Finally, the Consequence is immediate from Lemma 2, 
since we have J^(%-y L{t ^- ) ~Te T )y L ^ i g T )[F+p}(x,t)d9 -> 

in L^M, X) by Lemma 2, and ./^(Tr^r) - 
7^ T y i ( t)flT ))[F < ' T -F](a;,t)dd -> due to H 3 and Lemma 
1. 



3. Proofs 

Proof of Proposition 1. Due to i?2 and H3, we have 
e^ T (x,i) = 

dt !l[T T y L {t,T) - T Te 3> L(t , Te) ](i^ + P )d9 + h, (CI) 
where 

= Jo [(«(*) -v(t-r))T T y L{t , T) - 
(v(t) - v(t - 9r))T T0 y L ^ Te) ]d x {F^ + P )d9. 

Therefore, we obtain 

F Te y Ht , T e)} [F l ' T + P ]d9 + %*h. (C2) 

The second integral ii on the r. h. s. of Eq. (|C1[) vanishes 
if v is constant. If v is not constant, it tends to zero 
when r — > in i 1 ([0, in view of i?2 and H 3 Hi). 

Then, since ^ 7 (t) is bounded (by 1), \& 7 * Ii — > in 
L°°([0, T],X), hence in ^([O,! 1 ],!) since T is finite. 



APPENDIX D: PROBABILITY CURRENT 

1. The Fickian case 

The particles flux (probability current) J- e ' T (x, t) of the 
process x t ' T represents the average net number of walkers 
crossing point x at time t. Its dispersive and advective 
components are 

T l > T (x,t) ^^ T {x,t) + v(t)P £ ^{x,t), (Dl) 

where T^(x^t) denotes the contribution of dispersive 
jumps, and P l ' T (x,t) is the density of the process x e t ' T . 
Since each walker performs one dispersive jump per time 
step t, we have 

T i,r f +Q ° P LT (x-y,t)-P^(x + y,t) ^y^ 

(D2) 

where the function $(y/£) = ipi(z)dz represents the 
probability that dispersive jump length is larger than 
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y. Then, recalling that ipe is symmetric, we can rewrite 
j c ^ ) T (x, t) = T l jj T + {x, t) — Tp T _{x, t), where 



2 | +oc (DP e n(xTy,t)-(DP i n(x,t) it{ y )d ^ 

(D3) 

with D = £ 2 /2t. Hence, the contribution of dispersive 
jumps to the probability current is expressed through 
convolutions, whose kernel has a form £~ a ~ 1 JC(y / '£). 
Under some assumptions on /C, such mappings have 
a limit when I — > 0, that is a derivative of order 
a [13, [H, [H, [5?J • For the case considered here, a = 1, 
the lemma below shows that T£ converges in the 
hydrodynamic limit to —d x DP. The lemma applies 
to Eq. (|D3|) if DP has a uniformly bounded derivative 
with respect to x, provided that also P e,T — P — > 
more rapidly than £ In domains limited by reflecting 
boundaries, the dispersive flux needs corrections with 
respect to Eq. (|D2|) on the small scale, due to particles 
bouncing back at the walls [58]. Nevertheless, Eq. @ 
still holds for symmetric ipg with a finite second moment. 
Also, care must be taken when dealing with absorbing 
boundary conditions if ipi has a diverging second mo- 
ment dsj. 

Lemma 3. Let $ be a differentiablc function, integrable 
over R + , positive and decreasing. Then, for any inte- 
grable function g whose derivative is uniformly bounded, 



g(x + £y) - g(x) 
I 



<5>{y)dy 



1 dg 
2d^ 



Finally, for J{£) we have \J(£)\ < M\ /+~ ®{y)dy\, 
with M denoting the sup norm of the derivative of g. 
This term vanishes in the hydrodynamic limit, in view of 
the above choice of A(£) . 

2. The case of dispersion with immobile periods 

An explicit expression for the particles flux can be de- 
rived also for random walks with immobile periods. We 
will do it within scenario (51). We further assume that 
P^[ T and f l ' T converge according to hypothesis H 3 of 
Appendix [Cj on the basis of the probability for a tagged 
particle to cross x to the left/right during a small time 
interval. 

Particles that cross point x towards the right during 
time interval [t, t + dt] must be mobile, and have spent a 
time t' £ [0, t] in the mobile period. Moreover, they may 
i) or may not ii) have completed the single dispersive 
jump. In the former case i), they spent exactly a time r 
in that period, that began at point x — y — L(t, r) between 
instants t — r and t + dt — r, if the jump length is larger 
than y. Collecting all positive contributions y > gives 
for case i) the probability 



%yL(t,T)lf' T +r](x- y, t)$(hdydt 



<$>'{y)y 2 dy 



Case ii) can not occur if v(t) < 0. For positive v(t), it 
corresponds to particles that entered the mobile period 
between points x — v(t)dt — L(t,t') and x — L(t,t'), for 
all values of t' € [0, r], which yields the probability 



pointwise when £ — > 0. 

This proposition appears in (58|, within a slightly 
different context. Moreover, since $(y) = ipi(z)dz, 
we have $'(y) = —<pi(y). 

Proof. Let us denote A(£) a function, such that A(£) — > 
+oo when I — > 0, with £A(£) — * 0. For instance, we can 
assume A(£) — £~ a , with < a < 1. Then, we have 



g(x + £y) - g{x) 



which can be written as T{£) + J(£), with 

Jo ty 



f %>y L {t,t>)[f >T + r}( X) t)dfv(t)dt. 
Jo 

Upon dividing by dt, we recognize v(t)P^ T (x, t), accord- 
ing to Eq. ([5]). If v(t) > 0, crossings towards the left 
correspond to dispersive jumps and we only have 

T T y L{ttT) [f' T + r](x + y, t)${ y j)dydt 



and 



Hence, in view of Eq. (|12[) the probability current will be 
given by the sum of two terms 



~ P^{x -y,t)- P^ix + y, t) s { y )dy+v{t)P ^ t)> 



which tends to —d x DP m + v(t)P m as shown above, and 



g(x + £y) -g(x) 



y$(y)dy. 



Then, I{£) —* ^ / °° y<fr(y)dy, which integrating by 

f Q + °°yHy)dy = -kf + 



parts yields L + °° y<$>{y)dy = -\ <S>'(y)y 2 dy. 



e^(x - y, t)*(|)dy - / £ Lt (x + y, m{ V -)dy. 



This latter expression vanishes, according to the propo- 
sition 1 of Appendix [C] which proves Eq. (fT6|) . 
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